# 高斯分布的主成份分析
import argparse
from typing import Dict
import numpy as np
import matplotlib.pyplot as plt

class Chp01Sec05S1(object):
    def __init__(self):
        self.name = ''

    @staticmethod
    def startup(params:Dict = {}) -> None:
        print(f'高斯分布主成份分析 v0.0.1')
        plt.rcParams['figure.figsize'] = [16, 8]
        xC = np.array([2, 1])      # Center of data (mean)
        sig = np.array([2, 0.5])   # Principal axes
        theta = np.pi/3            # Rotate cloud by pi/3
        R = np.array([[np.cos(theta), -np.sin(theta)],     # Rotation matrix
                    [np.sin(theta), np.cos(theta)]])
        nPoints = 10000            # Create 10,000 points
        X = R @ np.diag(sig) @ np.random.randn(2,nPoints) + np.diag(xC) @ np.ones((2,nPoints))
        print(f'X: {X.shape};')
        fig = plt.figure()
        ax1 = fig.add_subplot(121)
        ax1.plot(X[0,:],X[1,:], '.', color='k')
        ax1.grid()
        plt.xlim((-6, 8))
        plt.ylim((-6,8))
        # 求PCA
        Xavg = np.mean(X,axis=1)                  # Compute mean
        B = X - np.tile(Xavg,(nPoints,1)).T       # Mean-subtracted data
        # Find principal components (SVD)
        U, S, VT = np.linalg.svd(B/np.sqrt(nPoints),full_matrices=0)
        print(f'B: {B.shape}; U: {U.shape}; S: {S.shape}; VT: {VT.shape}; Xavg: {Xavg.shape};')
        ax2 = fig.add_subplot(122)
        ax2.plot(X[0,:],X[1,:], '.', color='k')   # Plot data to overlay PCA
        ax2.grid()
        plt.xlim((-6, 8))
        plt.ylim((-6,8))
        #
        theta = 2 * np.pi * np.arange(0,1,0.01)
        # 1-std confidence interval
        v001 = np.array([np.cos(theta),np.sin(theta)])
        Xstd = U @ np.diag(S) @ np.array([np.cos(theta),np.sin(theta)])
        print(f'theta: {theta.shape}; v001: {v001.shape}; Xstd: {Xstd.shape};')
        # 绘制区间范围
        ax2.plot(Xavg[0] + Xstd[0,:], Xavg[1] + Xstd[1,:],'-',color='r',linewidth=3)
        ax2.plot(Xavg[0] + 2*Xstd[0,:], Xavg[1] + 2*Xstd[1,:],'-',color='r',linewidth=3)
        ax2.plot(Xavg[0] + 3*Xstd[0,:], Xavg[1] + 3*Xstd[1,:],'-',color='r',linewidth=3)
        # Plot principal components U[:,0]S[0] and U[:,1]S[1]
        ax2.plot(np.array([Xavg[0], Xavg[0]+U[0,0]*S[0]]),
                np.array([Xavg[1], Xavg[1]+U[1,0]*S[0]]),'-',color='cyan',linewidth=5)
        ax2.plot(np.array([Xavg[0], Xavg[0]+U[0,1]*S[1]]),
                np.array([Xavg[1], Xavg[1]+U[1,1]*S[1]]),'-',color='cyan',linewidth=5)
        plt.show()


def main(params:Dict = {}) -> None:
    Chp01Sec05S1.startup(params=params)

def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser()
    parser.add_argument(
        '--run_mode', action='store',
        type=int, default=1, dest='run_mode',
        help='run mode'
    )
    return parser.parse_args()

if '__main__' == __name__:
    args = parse_args()
    params = vars(args)
    main(params=params)